Screening in graphene antidot lattices 
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We compute the dynamical polarization function for a graphene antidot lattice in the random- 
phase approximation. The computed polarization functions display a much more complicated struc- 
ture than what is found for pristine graphene (even when evaluated beyond the Dirac-cone ap- 
proximation) ; this reflects the miniband structure and the associated van Hove singularities of the 
antidot lattice. The polarization functions depend on the azimuthal angle of the q-vector. We 
develop approximations to ease the numerical work, and critically evaluate the performance of the 
various schemes. We also compute the plasmon dispersion law, and find an approximate square-root 
dependence with a suppressed plasmon frequency as compared to doped graphene. The plasmon 
dispersion is nearly isotropic, and the developed approximation schemes agree well with the full 
calculation. 
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I. INTRODUCTION 

Graphene Antidot Lattices (GALs) have been pro- 
moted as a flexible platform for creating a tunable band 
gap thereby possibly allowing the realization of a number 
of technological applications^. In its basic form the GAL 
consists of a graphene sheet with a periodic nanometer 
scale perforation. Alternatively, GALs can be viewed as 
a special realization of a superlattice imposed on pris- 
tine graphene. Such a superlattice can be fabricated 
by a variety of technologies, for example by a periodic 
chemical modification, by selective adsorption of atoms 
or molecules on graphene, by periodic arrangement of 
electrostatic gates, or by an intrinsic or extrinsic regular 
corrugation. There is already a substantial literature on 
this general topic (see, e.g., Refs. 0t3|); here we focus on 
the antidot lattice case. Theoretically the properties of 
triangular GALs have been examined already quite sub- 
stantially (e.g., optical properties^, excitons^, electronic 
propertie d 10 ' 11 , electron-phonon couplin g 12 ' 13 , detection 
of edge states^, or details of band gap scaling^) . Most 
importantly, the experimental techniques for fabricating 
GALs have evolved rapidlyi£~— , presently reaching lat- 
tice constants down to a few tens of nanometers, where 
many of the interesting quantum mechanical effects pre- 
dicted by theory should become detectable. Of particu- 
lar interest are transport measurements, which we expect 
soon to be available. 

Most of the early works on GALs have focused 
on single-particle effects, such as the miniband struc- 
ture, sometimes supported by ab initio calculations 
using density-functional theory (DFT), or by some 
other less accurate but computationally more effective 
method, such as the density-functional based tight- 
binding (DFTB) 20 . However, it is necessary also to 
consider electron-electron interactions since these affect 
many physical properties, such as the charge carriers' 
interaction with light or with lattice vibrations, dielec- 
tric screening and plasmons, or the transport properties. 



The dynamical polarization function is a central object 
in these considerations, and it has been a topic of wide 
interest. In the literature expressions relevant for the 
polarization function of graphene, evaluated in the ran- 
dom phase approximation, can be found already from 
the "pre-graphene" era (see, e.g., Refs. [2 lll22l ]) . More 
recently, thorough studies of the polarization function 
in the Dirac-cone approximation have been reported in 
Refs. [23,24], but only very recently Stauber et al.— ^ 
gave similar results for full graphene dispersion at finite 
chemical potential. It should be noted that the accuracy 
of the RPA for pristine graphene is a highly nontrivial 
issue, see, e.g., Ref. [27|, and references cited therein. 
Many of the complications encountered in these stud- 
ies originate from the high-momentum cutoff needed in 
the Coulomb interaction, when extending the discussion 
beyond RPA, and are specific to the Dirac cone disper- 
sion law. Here, we consider the full GAL dispersion with 
a gap, and believe that the RPA is a justifiable starting 
point. Our calculations are conceptually straightforward: 
indeed the formal expressions for the Lindhard function 
for a general tight-binding model (which will be our un- 
derlying antidot lattice Hamiltonian) are standard text- 
book material (see, e.g., Sect. 8.5 in Ref. [28jj ) . Some 
complications inevitably arise because the antidot lattice 
unit cells contain tens or hundreds of atoms, in contrast 
to just two in pristine graphene. Analytic expressions 
can hardly be expected for such a complicated system, 
and much of our effort goes into examining whether sim- 
pler yet reasonably accurate models can reproduce the 
results of the full numerical calculations. The GALs can 
be realized with many different lattice symmetries, each 
of which may have some interesting physics of its own 
(see, e.g. Refs. 11.29.33]): whether a gap arises in the 
energy spectrum depends in a sensitive way on the sym- 
metry of the antidot lattice, and on the orientation of the 
antidot lattice with respect to the underlying graphene 
lattice. Here we examine one particular GAL (with a tri- 
angular lattice symmetry) in detail as a demonstration 
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of the theoretical methodology, and defer the discussion 
of other lattice symmetries until experiments performed 
on well-characterized samples become available. 

This paper is organized as follows. In Sect. [IT] we briefly 
introduce the basic models used in this work: the tight- 
binding representation of the antidot lattice states, and 
a simplified model, the "gapped graphene". In Sect. IIIII 
we give both analytical and numerical results for the po- 
larization function, address certain issues in the numeri- 
cal computation, and present a comparison between the 
various approximation schemes. Sect. IIVI addresses the 
plasmonic properties of GALs, and we end with a brief 
conclusion. 



II. THE MODELS 

A. Antidot lattice energy spectrum 

We use a tight-binding model to compute the band 
structure of the GALs. The wave function is written as 

M 

where (f>(r) is the normalized atomic wave function of the 
2p z orbital of carbon, and the j-sum runs over the atoms 
in the unit cell of the antidot lattice. R runs over the 
unit cells in the antidot lattice. The coefficients c^ k and 
the eigenenergies e n k are then obtained from the M x M 
eigenvalue problem, for each value of the wave vector k 
in the 1st Brillouin zone of the antidot lattice. Each 
GAL is characterized by a pair of indices {L, R}, where 
the L index gives the side length of the hexagonal unit 
cell, and the R index is the (approximate) radius of the 
hole, in units of the graphene lattice constant y/3a, where 
a = 1.42 A is the C-C distance^. Figure [1] shows an 
example of a typical GAL band structure. 



B. Gapped graphene 

The gapped graphene model is a phenomenological 
model used to describe band structures with a gap. 
Consider first the pristine graphene dispersion e° k = 
nf|</>(k)|, where t — 3.03 eV is the hopping integral, 
n = ±1 and 

<j)(k) = e lkSl + e^* 2 + e l]iS3 , (2) 

where Si are the three neighbors of the reference atom 
at the origin (<5i = a/2(-l,V3), S 2 = a/2(-l,-V3), 
S 3 = a(l,0)). Next, introduce a shift ±A on the onsite- 
energies for the carbon atoms on the A- and B- sublat- 
tices, respectively. This gives rise to a new spectrum, 

e„ k = nVA 2 +^ k | 2 , (3) 




FIG. 1: (Color online) Band structure of a {12, 3} GAL with 
a circular antidot with an armchair edge. The inset shows 
the atomistic configuration of the unit cell. The dashed line 
indicates the value of the chemical potential, fi = 0.05i (t is 
the hopping integral), used in most of the subsequent calcu- 
lations. 

where n = ±1, and corresponding eigenvectors 

( -nttj> k \ 
y/|e„ k |(| e , lk |-n A) \ , , 

)■ w 

Since </>(k) vanishes at the Dirac points K, K', the energy 
gap is Eg = e+.K — £-,k = 2|A|, and A can be chosen so 
that the spectrum fits the gap of a given antidot lattice 
spectrum. The gapped graphene can also be viewed as 
a minimal perturbation of the pristine graphene, e.g., in 
the case of very small antidots compared to the full unit 
cell. If, moreover, the chemical potential is close to the 
band edge so that the low-energy part of the band struc- 
ture dominates, we expect the gapped graphene model 
to be accurate. Below we give explicit examples of what 
"low-energy" means in practice. 



III. POLARIZATION FUNCTION 

A. Analytical results 

The polarization function is evaluated from the 
density-density response function for noninteracting elec- 
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trons: 



1 

nn'k 

»f(e n k) - ^F(en'k') 
e«k - e«'k' + w + irj ' 



dk|(nk|e- lq - r |n'k') 



(5) 



Here, tip is the Fermi-Dirac distribution, which we eval- 
uate at a fixed temperature of ksT = O.Olt through- 
out. In writing Eq. ([5]) we have assumed that local field 
effects can be neglected, i.e., the relevant wave vectors 
q are small compared to the reciprocal lattice vectors 
G oc 1/L. This allows us to work with a scalar dielectric 
function, rather than the more general object defined in 
the G, G'-space. The price is that certain collective in- 
tervalley modes will not be included (see, e.g., Ref. ;31J). 

In addition to the energy spectrum one also needs to 
evaluate the matrix element. In the following we list the 
required matrix element for the various cases of interest, 
including some well-known results for completeness. 



1. Pristine graphene 
The matrix element is2£ 



ikle-^ln'k') 



X— ( I - nn'ttc 



Jk',k+q 



y k+q 



(6) 



\<Pk\ |<Pk+q| 

where the various quantities are defined in Eq.@. 

2. Dirac cone 

With the spectrum linearized in the vicinity of the 
Dirac points, 



Sat 
n—k, 



(7) 



Eq.© simplifies, and one find; 



,23.24 



|(nk|e- iq >'k')| 2 =<5^ k+q ^[l + nn'cos(0 k , k+q )], (8) 
where # k , k + q is the angle between k and k + q. 



{L,R}={7,l},M/t = 1, 0=30" 
I <? |S = 0.01/(131,) 




0.3 0.6 0.9 1.2 1.5 

co/t (10" 3 ) 



FIG. 2: (Color online) Imxo(q, w) (top panel), and Rex5(q, w) 
(bottom panel) as a function oiuj/t for a {7, 1} GAL. <j> defines 
the azimuthal angle of the q- vector. The curve labeled "full" 
is evaluated using the full expression Eq. (|10p . while the curve 
labeled "diagonal" employs the approximation Eq. 



4- Antidot lattice 

Using the tight-binding wave functions ([T]) the matrix 
element becomes 



(nk|e^ 4q - r |n'k') 



3k' k 



M 



(4,)' 



n'k+q e 



(10) 

It is readily verified that the general expression (jlOl) re- 
produces the standard results (i.e., set M = 2, di = 0, 
and d2 = 83 for pristine graphene) . The challenge in the 
numerical applications is that for each q-point a large 
number of terms needs to be considered, and that a full 
Brillouin zone summation is required for every term. Be- 
low we discuss approximate methods of how to bypass 
this difficulty. 



B. Numerical results 



1. The limit qL < 1 



3. Gapped graphene 

Using the eigenenergies §5§ and eigenvectors ((U) one 
finds 

|(nk|e-^Vk')| 2 =<5 k ',k +q 
1 A 4 |0k| 2 |0k'| 2 + (kk| - nA) 2 (\e nk ,\- n'A) 2 



nn'Re 



2|e nk |(|e„ k | - nA)|e nk '|(|e nk '| - n'A) 



(9) 



As remarked above, in this work we do not consider 
local field effects, i.e., we restrict our discussion to wave- 
vectors that are small compared to the reciprocal lat- 
tice vectors, qL <C 1 (a similar restriction was used in 
Ref. [25| ) . We may exploit this restriction further by not- 
ing that the vectors dj entering the matrix element (flQ|) 
satisfy dj < L, and therefore the phase-factor can be ap- 
proximated as I exp(— iq ■ dj)\ ~ 1. Further, since the 
wave vectors q in general are small, for most of the k- 
vectors entering the Brillouin zone summations the eigen- 
vectors corresponding to k and k + q do not differ much, 
and therefore we can appeal to the orthonormality (in 
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{L,R)=(12,3}, p/t = 0.05 
| 9 |« = 0.01/(f3L) 




ffl/f (10 ) 



FIG. 3: (Color online) Imx5(q, w) (top panel) and Rex5(q, i^) 
(bottom panel) for the {12,3} GAL as a function of ui/t, at 
chemical potential \ijt — 0.05. The solid lines give the po- 
larization function (using the approximation pip ) for various 
azimuthal angles <j> while the dashed line gives the gapped 
graphene result. 



the band indices n, n) of the eigenvectors, and thereby 
obtain the final approximation for the matrix elements: 



(nk|e^ lq - r |n'k') 



?k\k- 



-q^n,n' 



fill 



We emphasize that since the limits of validity of the 
above expression are not rigorously given, one must care- 
fully check its reliability in each case. We have performed 
a large number of such numerical tests, and representa- 
tive results for a typical worst case scenario are given in 
Fig. [21 The results for the imaginary part, Fig. [5] (top 
panel) agree qualitatively well: all major features are re- 
produced. For the large limit, w/t > 1.25 in Fig. [2] 
a small systematic difference is observed (barely observ- 
able within the resolution of the figure); we have not been 
able to pinpoint the reason for this tiny difference. The 
results for the real part, Fig. [2] (bottom panel), were ob- 
tained by performing a Kramers-Kronig transformation 
of the corresponding imaginary part. Here we see per- 
haps an even better agreement between the full numerics 
and the approximate result. In general, the matrix ele- 
ment appears to be a smooth function of q, which does 
not introduce additional sharp features in the polariza- 
tion function. A particularly interesting observation is 
that the high-energy tails of the Rexo, which are impor- 
tant in the evaluation of the plasmon spectrum, are in 
near quantitative agreement, and we shall exploit this 
fact in Sect. HVl where we analyze the plasmon dispersion 
laws. 



2. Brillouin-zone summation 

The integrals over the Brillouin zone were carried out 
with an improved triangle method^, see the Appendix of 
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FIG. 4: (Color online) Imaginary part of the polarization 
function for the {12, 3} GAL as a function of uj/t, at chemical 
potential fi/t = 0.5. The solid lines give the full polarization 
function (using the approximation fTT}) for various azimuthal 
angles <f) while the dashed line gives the corresponding gapped 
graphene results. 



Ref. Q for details of our specific implementation. Very 
briefly, the method consists of generating a grid of k 
points, and writing the BZ-integrals as sums over tri- 
angles spanned by the grid points, and performing an 
analytic integration within each triangle using a linear 
approximation for the integrand. To evaluate the inte- 
grand, the M x Ai-eigenvalue and -vector problem must 
be solved for each k and k+q pair. By checking first that 
the Fermi-function difference in ([5]) is nonzero, as well 
as that the delta-function has a zero, a certain fraction 
of the eigenvalue problems could be eliminated. Fortu- 
nately, the method is rapidly converging as the k-mesh 
is made denser, as could be verified by bench-marking 
it against cases where analytic results are available, e.g., 
the Dirac cone dispersion. Typically we used ~ 5000 
k-points in the mesh. In our calculations we always eval- 
uate the imaginary part of Eq.([5]) first, and then compute 
the real part by a Kramers-Kronig transformation. 



3. An example: The {12, 3}-GAL 

Figure [3] gives the results for the imaginary and real 
part of the polarization function for the {12,3} GAL of 
Fig. 1, at a relatively low chemical potential, (i/t = 0.05. 
We have chosen to study this particular GAL in detail, 
and it is relevant to inquire how generic our results are. 
Of special interest is the effect of the shape and the edge 
of the antidot, for example whether the edge is of the 
armchair or of the zigzag type. It is well-known that 
non-dispersive zero-energy states may occur for certain 
structures (see, e.g., the examples shown in Refs. [TTIl2b] )]. 
which might lead to characteristic features in Xo(li u) )- 
Not much is known of the robustness of these edge states 
against disorder, or other edge reconstructions; see how- 
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ever a recent related study by Li et al22. A full study 
of these effects is beyond the scope of the present paper, 
where we want to illustrate the basic features of GALs 
with a relatively simple band structure, as the one illus- 
trated in Fig. 1. Another attractive feature of the chosen 
example is that the triangular GAL obeys a simple scal- 
ing law for the band gap 1 , which implies that our results 
are relevant also for larger structures, for which full nu- 
merical calculations would be extremely challenging. 

It is interesting to compare these results with those 
pertaining to pristine graphen o 25 ' 26 . A certain overall 
qualitative similarity persists. The imaginary part in- 
creases as a function of the frequency, and after reaching 
its maximum value it drops rapidly to zero. Superposed 
are certain sharp features: these are related to the van 
Hove singularities occurring at the band edges. Since 
there are many more bands in the GALs in a given en- 
ergy range than there are for pristine graphene, it is not 
surprising that the number of sharp peaks is larger for 
GALs. A similar qualitative resemblance can be seen for 
the real part of the polarization function: at low energies 
the real part plunges deeply (passing zero), whereafter it 
asymptotically approaches zero from the negative side. 
Again, the GAL exhibits much more fine structure. 

Also shown in Fig. [3] are the results for the gapped 
graphene model. We observe that the results depend 
on the azimuthal angle when the full GAL dispersion 
is used, whereas the gapped graphene model does now 
show such dependence. The overall shape of the gapped 
graphene model mimics fairly well the GAL results, but 
is obviously much smoother since the gapped graphene 
model with its only two bands cannot reflect the van Hove 
structure of the full GAL dispersion. 

While the gapped graphene model performs reasonably 
well at low chemical potentials, as expected, it fails even 
qualitatively at higher \x. An example is given in Fig.Q, 
where we show the imaginary part of the polarization 
function for fi/t = 0.5, and, as can be seen, there is 
only a token of resemblance between these results. The 
situation is even worse for the real part (not shown) . We 
conclude that if one is interested in the fine details of 
the polarization function, the gapped graphene model is 
reliable only at very low chemical potentials, and even 
there caution should be exercised. 



IV. PLASMONS 

We next examine the undamped plasmon modes for 
the {12, 3} GAL. As usual, the modes w(q) are identified 
as solutions to 



l-W(q)Rex5(q 1 w(q))=0, 



(12) 



where W(q) is the Coulomb interaction. We solve ([12")) 
numerically: for a fixed q we find the w p (q)'s which cause 
Eq. (IT2l to vanish, and use Imx(q, w p (q)) = as a crite- 
rion to identify the undamped modes. In all cases consid- 
ered the solution was unique. We shall compare the nu- 
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FIG. 5: (Color online) Plasmon dispersion for ti/t — 0.05. 
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FIG. 6: (Color online) Plasmon dispersion for fj,/t = 0.5. 



merically found GAL plasmons to the well-known disper- 
sions for pristine graphene 2 ^ w gr , and gapped graphene^ 
w g _ gr : 



W gr / 'fl = C 



w g-gr/^ 




(13) 



where C — y/ g s g v e 2 / (8ire r eoVF) with g s and g v spin- 
and valley degeneracies, respectively, and e r = 2.5, the 
relative dielectric constant for graphene on a Si02 sub- 
strate. In the general case the plasmon dispersion may 
depend on the direction of the wave vector q. However, 
as can be seen in Fig. [3] (bottom panel), at large values 
of ujjt - which determine the undamped plasmons be- 
cause Im%g vanishes there - the azimuthal dependence 
is very weak, and can be neglected. This result also fol- 
lows analytically from the definition of Xqi Eq.([5]): in 
the small qa limit it is possible to neglect the small en- 
ergy difference e„k — e„'k< as compared to uj, and the 
azimuth-independent l/cj-behavior of Re%Q follows, as 
seen in the numerics of Fig. 3, bottom panel. In Figs. [5] 
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and [6] we compare the plasmon dispersion laws for pris- 
tine graphene within the Dirac cone approximation, a 
gapped graphene with full 0(k), and the {12,3} GAL. 
Qualitatively, the GAL plasmon dispersion always lies 
below the other models. At low chemical potential the 
gapped graphene plasmon dispersion is essentially iden- 
tical with the full calculation (Fig. [5]), while for a higher 
/i (Fig. the gapped plasmon model works worse, as 
expected. 

V. CONCLUSIONS 

We have analyzed screening in graphene antidot lat- 
tices (GAL) within the random phase approximation. A 
general procedure for calculating the polarization func- 
tion is outlined. An efficient long wavelength approxima- 
tion is introduced, which eases the numerical task con- 
siderably, and the accuracy of this approximation is an- 
alyzed in terms of numerical examples. We also consider 
another phenomenological model: gapped graphene. We 
have chosen a {12,3} GAL as a generic system, and 
present results for this structure evaluated within the 
various approximation schemes. We conclude that the 
gapped graphene model works reasonably well for low 



chemical potentials (/i/t ~ 0.05), but that it fails even 
qualitatively for high chemical potentials {n/t ~ 0.5). 
We also determine the undamped plasmons for the 
{12,3} GAL. The plasmon dispersion law for the GAL 
has the same square-root behavior as pristine graphene, 
but it is significantly suppressed. The gapped graphene 
model works very well for low dopings, because the high- 
energy tails of the real part of the polarization are largely 
model independent. Future issues that need to be ad- 
dressed include the effects due to the GAL symmetry, 
and the geometry and edge structure of the nanoperfo- 
rations defining the antidots. 
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Note added in proof. During the technical process- 
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tailed account of the current-current correlation function 
of gapped graphene. Our results agree with theirs, when- 
ever there is overlap. 
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